import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import math
import random
import matplotlib as mpl
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn import metrics
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler,LabelEncoder
import plotly.graph_objs as go
import plotly.offline as offline
offline.init_notebook_mode()
import warnings
warnings.filterwarnings('ignore')
https://www.kaggle.com/crawford/principle-component-analysis-gene-expression/
Datos usados para clasificar pacientes con acute myeloid leukemia (AML) and acute lymphoblastic leukemia (ALL).
Golub et al "Molecular Classification of Cancer: Class Discovery and Class Prediction by Gene Expression Monitoring"
There are two datasets containing the initial (training, 38 samples) and independent (test, 34 samples) datasets used in the paper. These datasets contain measurements corresponding to ALL and AML samples from Bone Marrow and Peripheral Blood. Intensity values have been re-scaled such that overall intensities for each chip are equivalent.
# Load data
testfile = 'genes/data_set_ALL_AML_independent.csv'
trainfile = 'genes/data_set_ALL_AML_train.csv'
labels = 'genes/genes.actual.csv'
X_train = pd.read_csv(trainfile)
X_test = pd.read_csv(testfile)
y = pd.read_csv(labels)
# 1) Remove "call" columns from training a test
train_keepers = [col for col in X_train.columns if "call" not in col]
test_keepers = [col for col in X_test.columns if "call" not in col]
X_train = X_train[train_keepers]
X_test = X_test[test_keepers]
# 2) Transpose
X_train = X_train.T
X_test = X_test.T
X_train
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 7119 | 7120 | 7121 | 7122 | 7123 | 7124 | 7125 | 7126 | 7127 | 7128 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Gene Description | AFFX-BioB-5_at (endogenous control) | AFFX-BioB-M_at (endogenous control) | AFFX-BioB-3_at (endogenous control) | AFFX-BioC-5_at (endogenous control) | AFFX-BioC-3_at (endogenous control) | AFFX-BioDn-5_at (endogenous control) | AFFX-BioDn-3_at (endogenous control) | AFFX-CreX-5_at (endogenous control) | AFFX-CreX-3_at (endogenous control) | AFFX-BioB-5_st (endogenous control) | ... | Transcription factor Stat5b (stat5b) mRNA | Breast epithelial antigen BA46 mRNA | GB DEF = Calcium/calmodulin-dependent protein ... | TUBULIN ALPHA-4 CHAIN | CYP4B1 Cytochrome P450; subfamily IVB; polypep... | PTGER3 Prostaglandin E receptor 3 (subtype EP3... | HMG2 High-mobility group (nonhistone chromosom... | RB1 Retinoblastoma 1 (including osteosarcoma) | GB DEF = Glycophorin Sta (type A) exons 3 and ... | GB DEF = mRNA (clone 1A7) |
| Gene Accession Number | AFFX-BioB-5_at | AFFX-BioB-M_at | AFFX-BioB-3_at | AFFX-BioC-5_at | AFFX-BioC-3_at | AFFX-BioDn-5_at | AFFX-BioDn-3_at | AFFX-CreX-5_at | AFFX-CreX-3_at | AFFX-BioB-5_st | ... | U48730_at | U58516_at | U73738_at | X06956_at | X16699_at | X83863_at | Z17240_at | L49218_f_at | M71243_f_at | Z78285_f_at |
| 1 | -214 | -153 | -58 | 88 | -295 | -558 | 199 | -176 | 252 | 206 | ... | 185 | 511 | -125 | 389 | -37 | 793 | 329 | 36 | 191 | -37 |
| 2 | -139 | -73 | -1 | 283 | -264 | -400 | -330 | -168 | 101 | 74 | ... | 169 | 837 | -36 | 442 | -17 | 782 | 295 | 11 | 76 | -14 |
| 3 | -76 | -49 | -307 | 309 | -376 | -650 | 33 | -367 | 206 | -215 | ... | 315 | 1199 | 33 | 168 | 52 | 1138 | 777 | 41 | 228 | -41 |
| 4 | -135 | -114 | 265 | 12 | -419 | -585 | 158 | -253 | 49 | 31 | ... | 240 | 835 | 218 | 174 | -110 | 627 | 170 | -50 | 126 | -91 |
| 5 | -106 | -125 | -76 | 168 | -230 | -284 | 4 | -122 | 70 | 252 | ... | 156 | 649 | 57 | 504 | -26 | 250 | 314 | 14 | 56 | -25 |
| 6 | -138 | -85 | 215 | 71 | -272 | -558 | 67 | -186 | 87 | 193 | ... | 115 | 1221 | -76 | 172 | -74 | 645 | 341 | 26 | 193 | -53 |
| 7 | -72 | -144 | 238 | 55 | -399 | -551 | 131 | -179 | 126 | -20 | ... | 30 | 819 | -178 | 151 | -18 | 1140 | 482 | 10 | 369 | -42 |
| 8 | -413 | -260 | 7 | -2 | -541 | -790 | -275 | -463 | 70 | -169 | ... | 289 | 629 | -86 | 302 | 23 | 1799 | 446 | 59 | 781 | 20 |
| 9 | 5 | -127 | 106 | 268 | -210 | -535 | 0 | -174 | 24 | 506 | ... | 356 | 980 | 6 | 177 | -12 | 758 | 385 | 115 | 244 | -39 |
| 10 | -88 | -105 | 42 | 219 | -178 | -246 | 328 | -148 | 177 | 183 | ... | 42 | 986 | 26 | 101 | 21 | 570 | 359 | 9 | 171 | 7 |
| 11 | -165 | -155 | -71 | 82 | -163 | -430 | 100 | -109 | 56 | 350 | ... | 185 | 642 | 32 | 137 | -81 | 672 | 208 | 25 | 116 | -62 |
| 12 | -67 | -93 | 84 | 25 | -179 | -323 | -135 | -127 | -2 | -66 | ... | 48 | 224 | 60 | 194 | -10 | 291 | 41 | 8 | -2 | -80 |
| 13 | -92 | -119 | -31 | 173 | -233 | -227 | -49 | -62 | 13 | 230 | ... | 213 | 583 | 3 | 530 | -39 | 696 | 302 | 24 | 74 | -11 |
| 14 | -113 | -147 | -118 | 243 | -127 | -398 | -249 | -228 | -37 | 113 | ... | 267 | 440 | 52 | 229 | -4 | 431 | 269 | 8 | 163 | -22 |
| 15 | -107 | -72 | -126 | 149 | -205 | -284 | -166 | -185 | 1 | -23 | ... | 120 | 722 | 20 | 332 | -5 | 195 | 59 | 31 | 116 | -18 |
| 16 | -117 | -219 | -50 | 257 | -218 | -402 | 228 | -147 | 65 | 67 | ... | 79 | 631 | -26 | 455 | -62 | 736 | 445 | 42 | 246 | -43 |
| 17 | -476 | -213 | -18 | 301 | -403 | -394 | -42 | -144 | 98 | 173 | ... | 241 | 1215 | 127 | 255 | 50 | 1701 | 1109 | 61 | 526 | -83 |
| 18 | -81 | -150 | -119 | 78 | -152 | -340 | -36 | -141 | 96 | -55 | ... | 186 | 573 | -57 | 694 | -19 | 636 | 205 | 17 | 127 | -13 |
| 19 | -44 | -51 | 100 | 207 | -146 | -221 | 83 | -198 | 34 | -20 | ... | 318 | 397 | -48 | 1939 | -18 | 538 | 90 | -50 | 333 | -24 |
| 20 | 17 | -229 | 79 | 218 | -262 | -404 | 326 | -201 | 6 | 469 | ... | 225 | 1020 | -110 | 209 | -51 | 1435 | 255 | 53 | 545 | -16 |
| 21 | -144 | -199 | -157 | 132 | -151 | -347 | -118 | -24 | 126 | -201 | ... | 103 | 595 | -12 | 36 | 26 | 208 | 113 | -8 | 22 | -22 |
| 22 | -247 | -90 | -168 | -24 | -308 | -571 | -170 | -224 | 124 | -117 | ... | 158 | 402 | 57 | 253 | -52 | 1010 | 405 | 19 | 270 | -27 |
| 23 | -74 | -321 | -11 | -36 | -317 | -499 | -138 | -119 | 115 | -17 | ... | 129 | 1058 | 140 | 176 | -22 | 617 | 336 | 9 | 243 | 36 |
| 24 | -120 | -263 | -114 | 255 | -342 | -396 | -412 | -153 | 184 | -162 | ... | 176 | 725 | 13 | 249 | 1 | 646 | 391 | 81 | 203 | -94 |
| 25 | -81 | -150 | -85 | 316 | -418 | -461 | -66 | -184 | 164 | -5 | ... | 138 | 392 | 8 | 506 | 24 | 1034 | 69 | 24 | 807 | -41 |
| 26 | -112 | -233 | -78 | 54 | -244 | -275 | -479 | -108 | 136 | -86 | ... | 190 | 678 | 77 | 2527 | -36 | 838 | 313 | 21 | 145 | -19 |
| 27 | -273 | -327 | -76 | 81 | -439 | -616 | 419 | -251 | 165 | 350 | ... | 120 | 816 | 45 | 62 | -71 | 583 | 677 | -1 | 208 | 10 |
| 34 | -20 | -207 | -50 | 101 | -369 | -529 | 14 | -365 | 153 | 29 | ... | 260 | 1009 | -55 | 139 | -57 | 834 | 557 | -12 | 335 | -65 |
| 35 | 7 | -100 | -57 | 132 | -377 | -478 | -351 | -290 | 283 | 247 | ... | 93 | 336 | -45 | 170 | 12 | 752 | 295 | 28 | 1558 | -67 |
| 36 | -213 | -252 | 136 | 318 | -209 | -557 | 40 | -243 | 119 | -131 | ... | 234 | 1653 | 67 | 486 | -88 | 1293 | 342 | 26 | 246 | 23 |
| 37 | -25 | -20 | 124 | 325 | -396 | -464 | -221 | -390 | -1 | 358 | ... | 146 | 486 | -32 | 334 | 35 | 1733 | 304 | 12 | 3193 | -33 |
| 38 | -72 | -139 | -1 | 392 | -324 | -510 | -350 | -202 | 249 | 561 | ... | 103 | 1121 | 102 | 330 | -112 | 1567 | 627 | 21 | 2520 | 0 |
| 28 | -4 | -116 | -125 | 241 | -191 | -411 | -31 | -240 | 150 | 24 | ... | 173 | 755 | -23 | 573 | 42 | 987 | 279 | 22 | 662 | -46 |
| 29 | 15 | -114 | 2 | 193 | -51 | -155 | 29 | -105 | 42 | 524 | ... | 173 | 492 | 54 | 277 | -13 | 279 | 51 | 6 | 2484 | -2 |
| 30 | -318 | -192 | -95 | 312 | -139 | -344 | 324 | -237 | 105 | 167 | ... | 225 | 737 | 63 | 472 | 33 | 737 | 227 | -9 | 371 | -31 |
| 31 | -32 | -49 | 49 | 230 | -367 | -508 | -349 | -194 | 34 | -56 | ... | 36 | 592 | 57 | 215 | -22 | 588 | 361 | -26 | 133 | -32 |
| 32 | -124 | -79 | -37 | 330 | -188 | -423 | -31 | -223 | -82 | 176 | ... | 348 | 938 | -15 | 433 | -2 | 1170 | 284 | 39 | 298 | -3 |
| 33 | -135 | -186 | -70 | 337 | -407 | -566 | -141 | -315 | 206 | 321 | ... | 209 | 634 | -58 | 375 | -23 | 2315 | 250 | -12 | 790 | -10 |
40 rows × 7129 columns
# 3) Clean up the column names for training data
X_train.columns = X_train.iloc[1]
X_train = X_train.drop(["Gene Description", "Gene Accession Number"]).apply(pd.to_numeric)
# Clean up the column names for training data
X_test.columns = X_test.iloc[1]
X_test = X_test.drop(["Gene Description", "Gene Accession Number"]).apply(pd.to_numeric)
X_train.head()
| Gene Accession Number | AFFX-BioB-5_at | AFFX-BioB-M_at | AFFX-BioB-3_at | AFFX-BioC-5_at | AFFX-BioC-3_at | AFFX-BioDn-5_at | AFFX-BioDn-3_at | AFFX-CreX-5_at | AFFX-CreX-3_at | AFFX-BioB-5_st | ... | U48730_at | U58516_at | U73738_at | X06956_at | X16699_at | X83863_at | Z17240_at | L49218_f_at | M71243_f_at | Z78285_f_at |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | -214 | -153 | -58 | 88 | -295 | -558 | 199 | -176 | 252 | 206 | ... | 185 | 511 | -125 | 389 | -37 | 793 | 329 | 36 | 191 | -37 |
| 2 | -139 | -73 | -1 | 283 | -264 | -400 | -330 | -168 | 101 | 74 | ... | 169 | 837 | -36 | 442 | -17 | 782 | 295 | 11 | 76 | -14 |
| 3 | -76 | -49 | -307 | 309 | -376 | -650 | 33 | -367 | 206 | -215 | ... | 315 | 1199 | 33 | 168 | 52 | 1138 | 777 | 41 | 228 | -41 |
| 4 | -135 | -114 | 265 | 12 | -419 | -585 | 158 | -253 | 49 | 31 | ... | 240 | 835 | 218 | 174 | -110 | 627 | 170 | -50 | 126 | -91 |
| 5 | -106 | -125 | -76 | 168 | -230 | -284 | 4 | -122 | 70 | 252 | ... | 156 | 649 | 57 | 504 | -26 | 250 | 314 | 14 | 56 | -25 |
5 rows × 7129 columns
# 4) Split into train and test
X_train = X_train.reset_index(drop=True)
y_train = y[y.patient <= 38].reset_index(drop=True)
# Subet the rest for testing
X_test = X_test.reset_index(drop=True)
y_test = y[y.patient > 38].reset_index(drop=True)
Realiza un análisis exploratorio de los datos (correlaciones entre sí y con las clases, distribuciones,...). Usa las técnicas y gráficos que te parezcan más representativos.
En primer lugar conviene notar que nuestro conjunto de datos tiene un total de 7129 características, lo que convierte el problema en un problema de alta dimensión.
X_train.shape
(38, 7129)
En segundo lugar vamos a observar si existen missing values en el conjunto de datos.
(X_train.isnull().any()==False).all()
True
Como podemos observar, los datos están completos. A continuación codificamos cada clase de modo que {0: 'ALL', 1: 'AML'}.
label_enc = LabelEncoder().fit(y_train['cancer'].values)
y_train_enc = label_enc.transform(y_train['cancer'].values)
y_test_enc = label_enc.transform(y_test['cancer'].values)
Ahora vamos a seleccionar aquellas variables cuya correlación con la etiqueta 'cancer' sea superior en valor absoluto al 60%. Se representarán las matrices de correlación de la clase y los atributos seleccionados, así como los histogramas y los diagramas de clases correspondientes, de cada característica agrupado por clase.
# Get correlations of all train data
data = X_train.copy()
data['cancer'] = y_train_enc
corr = data.corr()
# Get columns with correlation greater than 60% in absolute value
corr_cancer = corr['cancer']
cond = abs(corr_cancer.values)>0.6
corr_cancer_filter = corr_cancer[cond]
# Create new data with these conditions
columns = list(corr_cancer_filter.index)[:-1]
data = X_train[columns]
data['cancer'] = y_train_enc
data.head()
| Gene Accession Number | AJ000480_at | D10495_at | D14874_at | D38128_at | D49950_at | HG1879-HT1919_at | HG3494-HT3688_at | HG4321-HT4591_at | J04027_at | L08177_at | ... | U41767_s_at | X58431_rna2_s_at | X83490_s_at | J03801_f_at | M19045_f_at | X14008_rna1_f_at | X16546_at | Z30644_at | U21689_at | cancer |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 895 | -208 | 163 | 36 | 75 | 400 | 223 | -64 | 96 | -110 | ... | 1057 | 424 | -66 | 921 | 382 | 557 | 174 | 1127 | -1660 | 0 |
| 1 | 749 | 334 | 176 | -952 | 129 | 528 | 1074 | -87 | 390 | 18 | ... | 785 | 368 | 30 | 4101 | 3606 | 2716 | 309 | 1378 | -973 | 0 |
| 2 | 1049 | -386 | 133 | -762 | 44 | 414 | 148 | -98 | 108 | -192 | ... | 1147 | 717 | 96 | 2799 | 2997 | 1716 | 226 | 1324 | -3091 | 0 |
| 3 | 1016 | 366 | 93 | -512 | 218 | 379 | -53 | -90 | 63 | -88 | ... | 1210 | 537 | 42 | 1166 | 1331 | 1143 | 151 | 1340 | -2041 | 0 |
| 4 | 634 | 408 | 75 | -897 | 110 | 324 | 353 | 174 | 357 | 30 | ... | 594 | 170 | 43 | 3250 | 3069 | 2917 | 672 | 806 | -976 | 0 |
5 rows × 85 columns
En segundo lugar se muestra la matriz de correlación de las 85 características seleccionadas junto con la variable cancer.
# Compute the correlation matrix
corr = data.corr()
# Generate a mask for the upper triangle
mask = np.triu(np.ones_like(corr, dtype=bool))
# Set up the matplotlib figure
f, ax = plt.subplots(figsize=(11, 9))
# Generate a custom diverging colormap
cmap = sns.diverging_palette(120, 630, as_cmap=True)
# Draw the heatmap with the mask and correct aspect ratio
sns.heatmap(corr, mask=mask, cmap=cmap, vmax=1, center=0,
square=True, linewidths=.5, cbar_kws={"shrink": .5})
plt.show()
Finalmente, se muestran los histogramas de cada característica agrupado por clases así como los diagramas de cajas.
fig = plt.figure(figsize = (20, 65))
j = 0
for i in data.columns[:-1]:
plt.subplot(21, 4, j+1)
j += 1
sns.distplot(data[i][y_train_enc==0], color='g', label = 'ALL')
sns.distplot(data[i][y_train_enc==1], color='r', label = 'AML')
plt.legend(loc='best')
fig.suptitle('Gene Data Analysis')
fig.tight_layout()
fig.subplots_adjust(top=0.95)
plt.show()
fig = plt.figure(figsize = (20, 65))
j = 0
for i in data.columns[:-1]:
plt.subplot(21, 4, j+1)
j += 1
ax = sns.boxplot(x='cancer', y=i, data=data)
fig.suptitle('Gene Data Analysis')
fig.tight_layout()
fig.subplots_adjust(top=0.95)
plt.show()
Observamos que en cada histograma las clases se solapan y no presentan una distribución conocida. También puede observarse la presencia de outliers en algunas características.
The analysis reveals that 21 principle components are needed to account for 80% of the variance. PC 1-3 add up to about ~33% and the rest is a slow burn where each component after PC8 contributes between 1-2% of the variance up until PC38 which is essentially zero. 1% is a decent amonut of variance and so the number of important PCs is up for interpretation.
# 5) Scale data
# (1) YOUR CODE HERE: Use the StandardScaler (separately for train and test sets)
scaler = StandardScaler().fit(X_train)
X_train_scl = scaler.fit_transform(X_train)
X_test_scl = scaler.fit_transform(X_test)
# 6) PCA Analysis and projection
components = 21
# YOUR CODE HERE:
# (2) Use PCA with this number of components on train set, with Y the result of the procedure
pca = PCA(n_components=components)
Y = pca.fit(X_train_scl)
# (3) Retrieve the explained variance ratio, and compute its accumulative sum
# save those values in variables var_exp and cum_var_exp
var_exp = Y.explained_variance_ratio_
cum_var_exp = np.cumsum(var_exp)
print(var_exp)
print(cum_var_exp)
[0.14987793 0.11977811 0.06600567 0.0488492 0.04632412 0.037219 0.03490875 0.03289664 0.02985074 0.0264505 0.02509445 0.02356326 0.0220301 0.02087031 0.01931271 0.01891243 0.01842475 0.01707946 0.01684773 0.01637575 0.01533421] [0.14987793 0.26965604 0.33566172 0.38451091 0.43083503 0.46805403 0.50296278 0.53585942 0.56571016 0.59216066 0.61725511 0.64081837 0.66284847 0.68371878 0.70303148 0.72194391 0.74036867 0.75744812 0.77429585 0.7906716 0.80600581]
Con 21 componentes se consigue un 80% de la varianza explicada.
Pregunta (1): ¿Qué pauta puede observarse en los valores de var_exp? ¿Cuál es la interpretación relativa de esos valores?
$\texttt{var_exp}$ es el vector que contiene la varianza total explicada por cada componente principal con respecto al conjunto. Vemos que cada vez la varianza explicada por cada componente principal es menor ya que las primeras componentes son las que maximizan la variabilidad de los datos.
Sea $\Sigma$ la matriz de covarianzas de los datos $X$ que queremos proyectar, la dirección de los componentes principales son los autovectores normalizados de $\Sigma$ correspondientes a los autovalores ordenador de mayor a menor. De manera que los autovalores $\lambda_1>\lambda_2>...$ de $\Sigma$ coinciden con la varianza de las proyecciones. Por tanto, las últimas componentes son aquellas cuyas proyecciones tienen menor varianza y eliminarlas no supone una gran pérdida de información.
# Plot the explained variance using var_exp and cum_var_exp
x = ["PC%s" %i for i in range(1,components)]
trace1 = go.Bar(
x=x,
y=list(var_exp),
name="Explained Variance")
trace2 = go.Scatter(
x=x,
y=cum_var_exp,
name="Cumulative Variance")
layout = go.Layout(
title='Explained variance',
xaxis=dict(title='Principle Components', tickmode='linear'))
data = [trace1, trace2]
fig = go.Figure(data=data, layout=layout)
offline.iplot(fig)
The first three components only explain 33% of the variance but we'll go ahead plot the projection to get a visual of it.
# Project first three components
Y_train_pca = pca.fit_transform(X_train_scl)
traces = []
for name in ['ALL', 'AML']:
trace = go.Scatter3d(
x=Y_train_pca[y_train.cancer == name, 0],
y=Y_train_pca[y_train.cancer == name, 1],
z=Y_train_pca[y_train.cancer == name, 2],
mode='markers',
name=name,
marker=go.Marker(size=10, line=go.Line(width=1), opacity=1))
traces.append(trace)
layout = go.Layout(
xaxis=dict(title='PC1'),
yaxis=dict(title='PC2'),
title="Projection of First Three Principle Components"
)
data = traces
fig = go.Figure(data=data, layout=layout)
offline.iplot(fig)
/home/maria/main/main/lib/python3.8/site-packages/plotly/graph_objs/_deprecations.py:378: DeprecationWarning: plotly.graph_objs.Line is deprecated. Please replace it with one of the following more specific types - plotly.graph_objs.scatter.Line - plotly.graph_objs.layout.shape.Line - etc. /home/maria/main/main/lib/python3.8/site-packages/plotly/graph_objs/_deprecations.py:434: DeprecationWarning: plotly.graph_objs.Marker is deprecated. Please replace it with one of the following more specific types - plotly.graph_objs.scatter.Marker - plotly.graph_objs.histogram.selected.Marker - etc.
Pregunta(2): Modificando la perspectiva de la figura con el ratón, ¿qué observas en cuanto a la separabilidad de las clases? Adjunta una imagen que apoye tus conclusiones.
En la imagen adjuntada, vemos que la separación de las clases es prácticamente lineal.

Realizar un análisis similar usando LDA, usando en este caso la información sobre el tipo de cancer de cada paciente. Puedes seguir la guía en https://www.apsl.net/blog/2017/07/18/using-linear-discriminant-analysis-lda-data-explore-step-step/
Es importante comentar que estamos trabajando con un problema de clasificación binaria y LDA utiliza una única componente ya que la dimensión del espacio vectorial en el que se proyectan los datos de entrenamiento es siempre menor o igual al número de clases. Por tanto, en la primera componente se explicará el 100% de la varianza del modelo y los datos se proyectarán en un vector normal al límite de decisión entre las dos clases.
sklearn_lda = LDA(n_components=1)
X_train_lda = sklearn_lda.fit_transform(X_train_scl, y_train_enc)
X_test_lda = sklearn_lda.transform(X_test_scl)
# Mean accuracy on the given training data and labels.
print(
"Mean accuracy on the given training data and labels: %s\n"
% str(sklearn_lda.score(X_train_scl, y_train_enc))
)
# Mean accuracy on the given test data and labels.
print(
"Mean accuracy on the given test data and labels: %s\n"
% str(sklearn_lda.score(X_test_scl, y_test_enc))
)
print(
"Explained variance ratio (first component): %s\n"
% str(sklearn_lda.explained_variance_ratio_)
)
for label in range(0,2):
mean = np.mean(X_train_lda[y_train_enc==label])
print('Mean Vector class %s: %s\n'
%(label, mean))
for label in range(0,2):
std = np.std(X_train_lda[y_train_enc==label])
print('Std Vector class %s: %s\n'
%(label, std))
Mean accuracy on the given training data and labels: 0.868421052631579 Mean accuracy on the given test data and labels: 0.6176470588235294 Explained variance ratio (first component): [1.] Mean Vector class 0: -0.5820367644856832 Mean Vector class 1: 1.4286356946466787 Std Vector class 0: 0.7736049956490063 Std Vector class 1: 1.3430444854571808
A continuación se realiza un histograma de las clases proyectadas de manera unidimensional. Para los datos de entrenamiento transformados vemos la clara separación entre ambas clases. Sin embargo, los datos de prueba se solapan y no se consigue un buen rendimiento, como se vió en la salida anterior del score.
def plot_hist(X, y, title):
sns.distplot(X[y==0], color='g', label = 'ALL')
sns.distplot(X[y==1], color='r', label = 'AML')
plt.legend()
plt.suptitle(title)
plt.tight_layout()
plt.xlabel('Projected feature')
plt.show()
plot_hist(X_train_lda, y_train_enc, 'Gene Data Train Analysis')
plot_hist(X_test_lda, y_test_enc, 'Gene Data Test Analysis')
Otro análisis similar se puede realizar proyectando sobre la dirección LD1 los datos transformados de entrenamiento y prueba con LDA. De nuevo se observa la separación entre clases para los datos de entrenamiento proyectados y el solapamiento para los datos de prueba proyectados. En la gráfica, también se ilustran las medias para cada etiqueta.
def plot_scikit_lda(X, y, title):
label_dict = {0:'ALL', 1:'AML'}
ax = plt.subplot(111)
for label, marker, color in zip(
range(0, 2), ('^', 's'), ('blue', 'red')):
X_train_lda_i = X[y == label]
X_train_lda_i_mean = np.array([X_train_lda_i.mean()])
plt.scatter(x=X_train_lda_i,
y=np.zeros(len(X_train_lda_i)), # flip the figure
marker=marker,
color=color,
alpha=0.5,
label=label_dict[label])
plt.scatter(x= X_train_lda_i_mean,
y=0, # flip the figure
marker='x',
color=color,
s=200)
plt.xlabel('LD1')
leg = plt.legend(loc='upper right', fancybox=True)
leg.get_frame().set_alpha(0.5)
plt.title(title)
# hide axis ticks
plt.tick_params(axis='both', which='both', bottom='off', top='off',
labelbottom='on', left='off', right='off', labelleft='on')
# remove axis spines
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.spines['left'].set_visible(False)
plt.grid()
plt.tight_layout
plt.show()
plot_scikit_lda(X_train_lda, y_train_enc, title='Train LDA via scikit-learn')
plot_scikit_lda(X_test_lda, y_test_enc, title='Test LDA via scikit-learn')
Utiliza k-means clustering con los datos originales y con los datos proyectados con PCA y LDA. ¿Qué observas?
Aunque K means sea un algoritmo no supervisado, este problema es binario, por tanto, vamos a forzar a que el número de clúster del algoritmo sea igual a 2. Así, podremos probar la eficiencia de Kmeans y ver si las predicciones de cada punto coinciden con la etiqueta real.
En primer lugar, se calcula la suma de los cuadrados de las distancias de los puntos desde sus respectivos centroides del clúster. Este valor, también llamado inertia, indica qué tan internamente coherentes son los clústeres
A continuación se realiza un análisis de siluetas. Este tipo de análisis sirve para estudiar la distancia de separación entre los clúster obtenidos y proporciona una forma de evaluar visualmente el número de clústers. Los coeficientes de silueta tienen un rango de $[-1,1]$ y los valores cercanos a $+1$ indican que la muestra está lejos de los clúster. Por su parte, un valor de 0 indica que la muestra está muy cerca de la frontera de decisión entre dos clúster. Finalmente, los valores negativos indican que esas muestras podrían haber sido asignadas al clúster equivocado.
# KMEANS
nclusters=2
kmeans = KMeans(n_clusters=nclusters, random_state=0)
X_train_km = kmeans.fit(X_train_scl)
y_train_km = kmeans.predict(X_train_scl)
y_test_km = kmeans.predict(X_test_scl)
print('Suma de cuadrados dentro del grupo: ', kmeans.inertia_)
# Silhouette analysis for data train
scores = metrics.silhouette_samples(X_train_scl, y_train_km)
sns.distplot(scores)
plt.title('Kmeans Data Train')
plt.show()
# Silhouette analysis for data test
scores = metrics.silhouette_samples(X_test_scl, y_test_km)
sns.distplot(scores)
plt.title('Kmeans Data Test')
plt.show()
Suma de cuadrados dentro del grupo: 241408.86104513207
Observando los resultados obtenidos vemos que tanto para los datos de prueba y entrenamiento, estos están muy cerca de la frontera de decisión de los dos clúster.
Como los datos tiene una gran dimensionalidad, vamos a realizar Kmeans sobre las dos componentes que maximizaban la varianza de los datos de entrenamiento. Para poder interpretar el rendimiento de Kmeans con PCA se mostrará tanto el valor de inertia así como los puntos clasificados por Kmeans respecto a las clases reales.
#PCA
pca = PCA(n_components=2).fit(X_train_scl)
X_train_pca = pca.transform(X_train_scl)
X_test_pca = pca.transform(X_test_scl)
#KMEANS
nclusters=2
kmeans = KMeans(n_clusters=nclusters, random_state=0)
X_train_pca_km = kmeans.fit(X_train_pca)
y_train_pca_km = kmeans.predict(X_train_pca)
y_test_pca_km = kmeans.predict(X_test_pca)
print('Suma de cuadrados dentro del grupo: ', kmeans.inertia_)
Suma de cuadrados dentro del grupo: 44273.60799143746
Vemos como ha disminuido el parámetro inertia de nuestro nuevo clasificador que utiliza los datos proyectados en las dos componentes principales. A continuación, se representan las componentes etiquedas por Kmeans para los datos de entrenamiento y prueba. También se ilustran las matrices de confusión para ambos datos.
def plotData(df1, groupby1, title1, df_2, groupby2, title2):
fig, (ax1, ax2) = plt.subplots(1, 2, figsize = (17,7))
# color map
cmap = mpl.cm.get_cmap('prism')
for i, cluster in df1.groupby(groupby1):
cluster.plot(ax = ax1, # need to pass this so all scatterplots are on same graph
kind = 'scatter',
x = 'PC1', y = 'PC2',
color = cmap(i/(nclusters-1)), # cmap maps a number to a color
label = "%s" % i,
s=30)
ax1.grid()
ax1.axvline(0, color='black')
ax1.set_title(title1);
for i, cluster in df_2.groupby(groupby2):
cluster.plot(ax = ax2, # need to pass this so all scatterplots are on same graph
kind = 'scatter',
x = 'PC1', y = 'PC2',
color = cmap(i/(nclusters-1)), # cmap maps a number to a color
label = "%s" % i,
s=30)
ax2.grid()
ax2.axvline(0, color='black')
ax2.set_title(title2);
# Train
X_pca = pd.DataFrame(X_train_pca, columns=['PC1','PC2'])
df_plot_train = X_pca.copy()
df_plot_train['ClusterKmeans'] = y_train_pca_km
df_plot_train['Cancer'] = y_train_enc
df_plot_train.head()
# Test
X_pca = pd.DataFrame(X_test_pca, columns=['PC1','PC2'])
df_plot_test = X_pca.copy()
df_plot_test['ClusterKmeans'] = y_test_pca_km
df_plot_test['Cancer'] = y_test_enc
plotData(df_plot_train, 'ClusterKmeans','Kmean cluster with PCA data train', df_plot_train, 'Cancer', 'Real labels with PCA data train')
plotData(df_plot_test, 'ClusterKmeans','Kmean cluster with PCA data test', df_plot_test, 'Cancer', 'Real labels with PCA data test')
print('Confusion matrix of train data:\n', metrics.confusion_matrix(df_plot_train['ClusterKmeans'],df_plot_train['Cancer']))
print('\nConfusion matrix of test data:\n', metrics.confusion_matrix(df_plot_test['ClusterKmeans'],df_plot_test['Cancer']))
Confusion matrix of train data: [[19 1] [ 8 10]] Confusion matrix of test data: [[13 9] [ 7 5]]
Tanto en los gráficos como en la matriz de confusión se ve claramente los verdaderos positivos y negativos, así como los errores de clasiciación, siendo mucho mayores para los datos de prueba.
De nuevo, se calcula el parámetro de inertia, se represenan los histogramas de los clúster proyectados en la componente de LDA así como la recta con los valores proyectados, y se calculan las matrices de confusión de las predicciones respecto a los valores reales.
## Kmeans with LDA
nclusters=2
kmeans = KMeans(n_clusters=nclusters, random_state=0)
X_train_lda_km = kmeans.fit(X_train_lda)
y_train_lda_km = kmeans.predict(X_train_lda)
y_test_lda_km = kmeans.predict(X_test_lda)
print('Suma de cuadrados dentro del grupo: ', kmeans.inertia_)
plot_hist(X_train_lda, y_train_lda_km, 'LDA Kmeans Data Train Analysis')
plot_hist(X_test_lda, y_test_lda_km, 'LDA Kmeans Data Train Analysis')
plot_scikit_lda(X_train_lda, y_train_lda_km, title='Train LDA Kmeans via scikit-learn')
plot_scikit_lda(X_test_lda, y_test_lda_km, title='Test LDA Kmeans via scikit-learn')
print('Confusion matrix of train data:\n', metrics.confusion_matrix(y_train_lda_km,y_train_enc))
print('\nConfusion matrix of test data:\n', metrics.confusion_matrix(y_test_lda_km,y_test_enc))
Suma de cuadrados dentro del grupo: 21.357177393370186
Confusion matrix of train data: [[26 4] [ 1 7]] Confusion matrix of test data: [[18 10] [ 2 4]]
En este caso, el parámetro de inertia es muy bajo, obteniendose un valor de $21.35$. Los histogramas representan claramente la separación entre clases, aunque recordemos que para el caso de los datos de prueba proyectados sobre la componente de LDA las clases continuaban solapadas. Finalmente, en la matriz de confusión se resumen los datos bien y mal clasificados de cada clase, habiéndose reducido considerablemente el número de falsos negativos para los datos de entrenamiento.